Island, pit and groove formation in strained heteroepitcLxy 
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We study the morphological evolution of strained heteroepitaxial films using a kinetic Monte 
Carlo method in three dimensions. The elastic part of the problem uses a Green's function method. 
Isolated islands are observed under deposition conditions for deposition rates slow compared with 
intrinsic surface roughening rates. They are semi-spherical and truncated conical for high and low 
temperature cases respectively. Annealing of films at high temperature leads to the formation of 
closely packed islands consistent with an instability theory. At low temperature, pits form via a 
layer-by-layer nucleation mechanism and subsequently develop into grooves. 

PACS numbers: 68.65.-k, 68.65.Hb, 81.16.Dn, 81.16.Rf 



Epitaxial growth techniques have been used to deposit 
strained coherent films on substrates of a different mate- 
rials with a mismatched lattice constant. This is called 
heteroepitaxy. Many experiments have shown that be- 
yond a threshold film thickness, an array of three dimen- 
sional (3D) nanosized islands self-assembles under favor- 
able growth conditions P, 0, 01- These results are of 
considerable interest since the islands behave as quan- 
tum dots and are expected to find applications in future 
microelectronic devices. The most intensively studied ex- 
amples include Ge/Si(100)and more generally its alloy 
variant Sii_^Ge^/Si(100) The island mor- 

phology depends strongly and often non-trivially on the 
lattice misfit dictated by the Ge concentration as well 
as growth conditions including temperature and deposi- 
tion rate. In addition, other interesting nanostructures 
including 3D pits, grooves and quantum dot molecules 
composed of coupled islands and pits are also generated 
under appropriate conditions [3, [l3 • 

In this letter, we report large scale 3D kinetic Monte 
Carlo simulations on the morphological evolution of 
strained layers. Our simulations generate morphologies 
very reminiscent of those observed under various growth 
or annealing conditions. We should note that the simula- 
tion of strained layers is computationally challenging due 
to the long range nature of elastic interactions. Previous 
atomistic simulations are limited to two dimensions (2D) 
[Til IT^ ITsl IT^ or sub- monolayer coverage 0. Con- 
tinuum computations are less difficult but cannot reli- 
ably account for faceted surfaces and fluctuations which 
are especially important at the early stage of roughening 

EliElil. 

We model the film and substrat e sy stem by a simple 
cubic lattice of balls and springs The sub- 

strate consists of 64 X 64 X 64 atoms. Periodic boundary 
conditions in lateral directions and fixed boundary con- 
ditions for the bottom layer are assumed. The substrate 
has a lattice constant as — 2. 72 A which gives an atomic 
density appropriate for crystalline silicon. The lattice 
constant a/ of the film is related to the lattice misfit 
e= {uf — as)/af. Nearest neighboring (NN) and next 



nearest neighboring (NNN) atoms are directly connected 
by elastic springs with force constants fci — 2eV/a'j. and 
^2 — ki respectively. The elastic couplings of adatoms 
with the rest of the system are weak and are completely 
neglected. 

Our algorithm imposes solid-on-solid conditions with 
atomic steps limited to at most one atom high. Every 
topmost atom in the film can hop to a different random 
topmost site within a neighborhood of I x I columns with 
equal probability. We put I = 33. Decreasing the hop- 
ping range does not alter our results significantly. The 
hopping rate r„i of a topmost atom m follows an Arrhe- 
nius form 



Tm = Rq exp 



»lm7l + "-2m72 - ^E„i - Ep 

knT 



(1) 



Here, rii^ and 7i2m are the number of NN and NNN of 
atom m respectively while 71 — 0.085ey and 72 = 7i/2 
are the corresponding bond strengths. The elastic energy 
of the hopping atom is denoted by AE„i and will be 
explained later. Finally, we put Eq = 0.415eV and Rq = 
2Dpl{aasf with Dq = 4.1 x lO^^A^-i and cr^ = ^2/6- 
This gives the a ppr opriate adatom diffusion coefficient 
for silicon (100) ^M,- Our choice of the ratios ki/k2 and 
7i /72 maximizes the isotropy of the system. 

The elastic energy, AEm, has to be repeatedly calcu- 
lated during a simulation; this dominates the CPU time. 
AE„i is defined as the difference in the strain energy Es 
of the whole lattice at mechanical equilibrium when the 
site is occupied minus that when it is unoccupied. Cal- 
culating Eg requires solving a long-range elasticity prob- 
lem to obtain the atomic positions of every atom in the 
film and the substrate. We have found it possible to 
significantly speed up the calculation by appling an ex- 
act Green's function method. A method of this type was 
introduced by Tewary (2Ql | in the context of point impuri- 
ties. We generalized the technique to free surfaces in Ref. 

. The result of these developments is that we can solve 
the elastic problem at a surface site using reduced equa- 
tions involving only other surface atoms. Moreover, we 
use a surface coarsening scheme in which morphological 




FIG. 1: Surface from simulation of deposition at lOOOK and 
20000 MLs~^ in top view (a) and 3D view (b). The gray scale 
shows the local height of the surface and the exposed part of 
the substrate is shaded in brown. 



details of the surface far away from atom m are averaged 
[T3I I. As a result, calculating AEm involves only about 
160 effective particles and takes less than one second on 
a 3GHz Pentium computer. Hopping events are then 
sampled using an acceptance-rejection algorithm aided 
by quick estimates of AEm which enables a high accep- 
tance probability. A simulation reported here typically 
involves 10® successful hopping events and takes 10 days 
to complete. We have considered large misfit and in some 
cases also high deposition rate so that the computations 
can be manageable. 

We have simulated deposition of films with 8% lattice 
misfit at temperature lOOOK and deposition rate 20000 
MLs~^. Figurenshows the resulting morphology from a 
typical run at a nominal film thickness of 3MLs. Isolated 
semi-spherical islands are observed. Most of them nucle- 
ate when the nominal coverage is about 1 ML and then 
grow steadily as more atoms are deposited. Coarsening 
via exchange of atoms among islands (Ostwald ripening) 
also occurs. Some small islands shrink and vanish even- 
tually. However, coalescence of islands is suppressed by 
their mutual elastic repulsion ,2lJ. In fact, the edges 
of neighboring islands are often deformed to avoid each 
others. 

In our simulations, as in experiment, the deposition 
rate has a substantial effect on surface morphology. At 
the rate considered above, island growth is limited by the 
supply of atoms. Individual islands have already relaxed 
to their equilibrium shapes. That is, deposition is slow 
relative to the formation dynamics and geometrical relax- 
ation of islands. In contrast, at lower deposition rates, 
we observe that islands become larger and less dense be- 
cause there is more time for coarsening. For deposition 
faster than island formation, layers of ator ns q uickly ac- 
cumulate before the resulting film roughens |l3| . With an 
abundant supply of atoms, we observe that the roughen- 
ing dynamics is similar to that for annealing except for a 
trivial vertical drift of the whole surface. We will discuss 
annealing next. 

We have simulated annealing of initially flat films with 
lOMLs of atoms and 6% lattice misfit at lOOOK. Fig- 
ures|2Ia)-(c) show snapshots of the evolution. 2D islands 
and pits first develop leading to a high step density [Fig. 
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FIG. 2: Snapshots from annealing of an initially flat film at 
lOOOK at time t = 20 (a), 50 (b) and 100/iS (c), and a plot of 
surface width w against t for 5 independent runs (d). 



|2ta)]. At this point, the film is still relatively fiat and 
highly stressed. The misfit has little impact on the mor- 
phology except for an enhancement of the step density 
due to a reduction of the effective step free energy. As the 
roughness increases, long-range elastic interactions begin 
to dominate and lead to the formation of 3D islands and 
pits with gentle slopes [Fig. |2Ib)]. Subsequently, well 
developed 3D islands bounded by a network of grooves 
emerge [Fig. |2Ic)]. Note that the surface inclination at 
many grooves has reached its maximum value allowed in 
our model. In experiments on the Si-Ge system, grooves 
are often bounded by [115] facets. The physical reasons 
for this might be the same, though our model is too crude 
to select among facets. 

The (100) surface studied above does not act like a true 
facet as is evident from the abundance of surface steps in 
Fig. Ola). This indicates that lOOOK is above the surface 
roughening transition temperature. Thus, the surface en- 
ergy varies smoothly with the local inclination. In this 
situation, the strain-induced roughening of an unfaceted 
surface should be described by the Asaro-Tiller-Grinfeld 
instability theory j23 | which predicts that random per- 
turbations of the surface at sufficiently long wavelength 
spontaneously amplify. The surface will gradually be 
dominated by modulations at the most unstable wave- 
lengths. 

Our annealing results at lOOOK are consistent with the 
instability theory. This is supported by a few characteris- 
tic features. First, the sidewalls of the newly emerging is- 
lands are gentle and their inclinations increase gradually 
rather than abruptly. Moreover, the island base areas 
stay relatively constant. As a further evidence. Fig. |2Id) 
plots the r.m.s. surface width w against the annealing 
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FIG. 3: Surface from simulation of deposition at 600K and 
10 MLs-\ 
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FIG. 4: Snapshots from annealing of an initially flat film at 
600K at time t = 0.1 (a), 0.15 (b) and 0.22 s (c), and a plot 
of w against t for 5 independent runs (d). 



time t for 5 independent runs. We observe that w in- 
creases steadily and the ensemble fluctuations are small 
as expected for barrierless processes. The morphologi- 
cal development also qualitatively resembles the initial 
evolution of Sii_a;Gc2;/Si(100) films at high temperature 
and low misfit (31 • Tersoff, et. al have argued that the 
Sii_a;Gea;(100) surface under these conditions is not a 
true facet and theories based on unfaceted surfaces 
should apply. 

Next, we consider a lower temperature, 600K, which 
gives drastically different morphologies indicating dis- 
tinct roughening mechanisms. Figure |2| shows a surface 
at a nominal coverage of 2MLs from a simulation of de- 
position at 8% misfit and 10 ML s~^. We again observe 
isolated islands but they now take the shapes of trun- 
cated cones. Most islands are out of equilibrium as their 
heights are clearly limited by significant energy barriers 
for upper layer nucleation. 

We have also simulated annealing at 600K. Figures 
131 a)- (c) show three snapshots from a typical run. A 
large 2D island and a few smaller 2D pits first appear 



[Fig. Efa)]. Later, 3D pits develop [Fig. Efb)]. They 
then become increasingly eccentric and gradually turn 
into grooves [Fig. ^c)]. Analogous 3D structures are 
also observed for deposition at rates fast compared to 
roughening. 

It is interesting to note that only part of the surface 
roughens even after a long annealing time in sharp con- 
trast to the high temperature case. This strongly sug- 
gests that the surface is a true facet and 600K is below 
the roughening temperature. The surface energy should 
be a singular function of the slope and instability theory 
should not apply. For this situation, a nucleation the- 
ory has been suggested for 3D island or pit formation on 
faceted surfaces 24, 25]. According to this approach, an 
island or pit has to overcome an energy barrier associated 
with a critical volume before it can be stable. Experi- 
ments on island formation at low temperature and high 
misfit have indicated better agreement with nucleation 
theory ^24] . 

Figure Sl^d) plots the r.m.s. surface width w against 
time from 5 independent runs during early stage of rough- 
ening. There are large ensemble fiuctuations supporting 
the relevance of nucleation processes. However, there ex- 
ists no dominating jump in w associated with a single 
successful nucleation event after which w grows steadily. 
Instead, multiple relatively rapid increments can be ob- 
served and are associated with the creation of lower lay- 
ers in the dominant pits. Conventional nucleation theory 
does not successfully describe the effects of the large bar- 
riers for nucleation of further layers. The formation of 
3D pits in Figs. Ela)-(c) and in fact also of the 3D is- 
lands in Fig. 13 should best be described by a sequence 
of layer-by-layer nucleation events. For a growing pit 
for instance, atoms are ejected continuously while lateral 
expansion takes place at constant pit depth. Once the 
bottom becomes sufficiently large, nucleation of a further 
layer will be possible. The growth is thus based on the 
correlated processes of continuous lateral expansion and 
periodic sudden nucleation of deeper layers. The associ- 
ated rates depend not only on the pit geometry but also 
on the presence of nearby islands or pits due to both ex- 
change of atoms and elastic interactions. Recent theories 
on the elasticity of step mounds should be particularly 
relevant for further analysis t26i.i27i|. 

The selection mechanism between islands and pits also 
deserves further explanation. Continuum elasticity the- 
ory shows that islands and pits with infinitesimal slopes 
relieve elastic energy equally well [24]. It is visually ap- 
parent that an up-down symmetry exists for surfaces in 
Fig. [2Ia) and to a lesser extent also in Fig. |2tb). How- 
ever, pits are increasingly favored energetically compared 
to islands as local slopes become steeper At low 

temperatures, the energy difference is already significant 
for single layered structures. Specifically, asymmetry be- 
tween 2D islands and pits is already apparent in Fig. 
^a). There is typically one dominant island but a few 
smaller pits. This is because the lower energy of pits 
also implies a lower nucleation barrier. Pits can hence 
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nucleate more quickly and are more abundant. Because 
new islands are not nucleated, the existing one absorbs 
all ejected atoms and grows quickly. Furthermore, the 
better stability of pits also explains the development of 
3D pits rather than 3D islands in Figs. Elb)-(c). We 
have observed 3D islands only in Fig. 13 for slow deposi- 
tion. This is because under this slow deposition rate, 3D 
islands are already able to develop before a thick enough 
film can be formed to accommodate the pits Ex- 
perimentally, the selected structure also turns from 3D 
islands to 3D pits upon lowering the temperature and 
increasing the deposition rate |^. 

An interesting transition from pits to grooves is also 
observed in Figs. Elb)-(c). For shallow pits, a square 
base is energetically preferred to a rectangular one 
This explains the more rounded shapes of the pits in Fig. 
01b). As the pits enlarge, their sidewalls also become 
steeper to relieve the stress more efficiently. Grooves are 
then energetically preferred to rounded pits because their 
linear extents are larger and can lead to stress relief over 
a much wider region. Formation of grooves in Fig. EJc) 
is further enhanced by the stress around a 2D island. 
This is closely related to the phenomena of cooperative 
nucleation |5( and also trench formation around 3D is- 
lands [2^ . The presence of a neighboring island however 
is not essential as we have also observed pits turning into 
grooves far away from any islands. Grooves are also ob- 
served in experiments from annealing of pits |lO| . 



In conclusion, we have applied a kinetic Monte Carlo 
method in 3D to study morphological structures gen- 
erated from deposition and annealing of strained het- 
eroepitaxy. Under deposition conditions, morphologies 
depend dramatically on whether deposition is slow com- 
pared to the intrinsic roughening rate of the surface as 
in the 2D case jl^ . For slow deposition, isolated islands 
result and their formation and development are limited 
by the supply of atoms. In contrast during fast deposi- 
tion, 3D structures form only after layers of atoms have 
accumulated and are similar to those from annealing of 
initially flat films. Morphologies from annealing further 
show strong dependence on temperature which deter- 
mines whether the initial surface is faceted. Upon an- 
nealing at high temperature, unfaceted surfaces develop 
arrays of 3D islands via the Asaro-Tiller-Grinfield insta- 
bility. In contrast, faceted surfaces at low temperature 
develop 3D pits via a layer-by-layer nucleation mecha- 
nism. The pits later turn into grooves. We suggest that 
the selection mechanisms between islands and pits as well 
as between pits and grooves are of energetic origin. Many 
of the general trends that we observe in our simulations 
are similar to experimental results. 
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